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Abstract 

The present work is devoted to the simulation of a strongly magnetized plasma 
considered as a mixture of an ion fluid and an electron fluid. For the sake of 
simplicity, we assume that the model is isothermal and described by Euler equations 
coupled with a term representing the Lorentz force. Moreover we assume that both 
Euler systems are coupled through a quasi-neutrality constraint of the form rij = n e . 
The numerical method which is described in the present document is based on 
an Asymptotic-Preserving semi-discretization in time of a variant of this two-fluid 
Euler-Lorentz model with a small perturbation of the quasi-neutrality constraint. 
Firstly, we present the two-fluid model and the motivations for introducing a small 
perturbation into the quasi-neutrality equation, then we describe the time semi- 
discretization of the perturbed model and a fully-discrete finite volume scheme based 
on it. Finally, we present some numerical results which have been obtained with 
this method. 
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1 Introduction. 

This paper is devoted to the construction of a numerical scheme for the simulation of a 
two-fluid Euler-Lorentz model: such a model represents the evolution of a mixture of an 
ion gas and an electron gas which are submitted to the Lorentz force. More precisely, we 
focus in this paper on a situation involving a strong Lorentz force and a low Mach number 
regime for both ion and electron fluids, i.e. we assume that pressure and Lorentz forces 
are of the same order as r _1 where r > is the square of the ion Mach number and also 
represents the ratio between the ion gyro-period and the characteristic time scale of the 
experiment. When r converges to 0, we reach an asymptotic regime which is referred to 
as the drift-fluid regime or gyro- fluid regime: in this limit regime, the pressure force for 
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the ions and for the electrons is balanced by Lorentz force. In return, the momentum 
equations within the two-fluid Euler system degenerate into a pair of equations in which 
the parallel components of the ion and electron velocities can be viewed as Lagrange mul- 
tipliers of the zero total force equations in the direction of the magnetic field (see [5l IT6]). 

Such a model describes plasma physics experiments involving strong external magnetic 
fields, such as Magnetic Confinement Fusion (MCF) experiments in tokamak reactors. In 
such a case, the rescaled gyro-period of the confined particles, which is denoted by r, can 
be close to 0. The assumption that r is very small leads to a singularly perturbed Euler- 
Lorentz model. The limit r — > is referred to as the gyro-fluid limit (see [U EH [30] ) . 
It is also possible to consider a gyro-kinetic approach when a kinetic model is considered 
from the onset, instead of fluid equations (see [31 SI EQl E21 EE1 EHl EH EH EHl E]). 
For generalities on asymptotic regimes for fusion plasmas physics, we refer to [13]. 

In many experimental cases, the value of r is not uniform and can vary a lot between 
subdomains of the tokamak. Additionally, it may depend on the time variable: in most 
of MCF experiments, r is very small in the plasma core whereas it can be of order 1 far 
from the plasma core. From a numerical point of view, the usual approach for simulating 
both cases together consists in a domain decomposition according to the local value of r. 
More precisely, we choose to simulate the initial r-dependent model in the regions where 
r = 0(1), and we choose the limit model in the regions where r< 1. Such an approach 
involves different numerical methods for solving either the Euler-Lorentz model or its 
drift-fluid limit according to the value of r. Generally, the coupling of these methods is 
not straightforward and presents several drawbacks such as the treatment of the interface 
position (or cross-talk region): indeed, it can depend on the time variable and, in most 
cases, costly algorithms are required to simulate the motion of the interface and to couple 
it with the space mesh. 

We choose a different method based on the resolution of the r-dependent Euler-Lorentz 
model and on the design of a scheme which is able to handle both the cases r = 0(1) 
and r < 1. Then this so-called Asymptotic-Preserving (AP) method provides consistent 
approximation of the Euler-Lorentz model when r = 0(1) and of its limit regime when 
r —¥ 0, and does not require a r-dependent stability condition. As a consequence, such a 
method can be used on the whole simulation domain for both the r = 0(1) and r < 1 
regimes. AP schemes have been introduced by S. Jin [32] and were applied on tranport 
models and their diffusive limits. Other applications can be found in plasma physics (see 
[21 [9j [TOj [TTJ [HJ [151 El H] f° r quasi-neutrality regimes and [5j [16] for strong magnetic 
fields regimes), low Mach number fluid dynamics (see [13 H2]), or other types of transport 
problems (see PElIHl[I21[211[251ESlEn]) or diffusion problems (see [T3]). 

The present paper has two main goals: the first one is to propose a new equivalent 
formulation of the two-fluid Euler-Lorentz model when r > 0. In this new formulation, 
the parallel velocity equation for r = explicitly appears as the limit of the parallel 
velocity equation for r > by contrast to the original formulation. This reformulation is 
the building block for the Asymptotic-Preserving scheme for the two-fluid Euler-Lorentz 
model. For the scheme being simple, we choose an isothermal pressure law for both ion 
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and electron fluids. We also assume that the fluids are coupled through a quasi-neutrality 
constraint which allows us to compute the self-consistent electric field within the Lorentz 
term, the magnetic field being external and given. This work is the last development 
of a program started in [5l [16]: in [5] and [16], the one- fluid isentropic Euler-Lorentz 
model with given electric and magnetic fields has been investigated (in [16] , under a uni- 
form magnetic field and in [5j, with any magnetic field and arbitrary coordinate system). 
Here, the specificity of this work is to consider a two-fluid Euler-Lorentz model with self- 
consistent electric field, computed through the quasi-neutrality hypothesis. This leads 
to a system of two coupled anisotropic diffusion equations, which brings some specific 
difficulties. In particular, it involves some singularity which will be treated through a reg- 
ularization procedure. This work also bears relations with [T3| [15] which are concerned 
with more general anisotropic diffusion equations (but not in the context of the Euler- 
Lorentz model). Again, [13] deals with a uniform anisotropy direction and [T5] with an 
arbitrary anisotropy direction and arbitrary coordinate systems compared to the direction 
of the anisotropy. 

The present paper is organized as follows. In section 2, we present the isothermal 
two-fluid Euler-Lorentz model and the drift-fluid limit regime. The section 3 is devoted 
to the reformulation of the r-dependent Euler-Lorentz model leading to a new equivalent 
formulation of these equations when r > which is equivalent to the drift-fluid limit 
when r = 0. In section 4, we present a time semi-discretization of the Euler-Lorentz 
model which is also consistent with the new formulation of the model. Since this time 
semi-discrete scheme involves an ill-posed diffusion problem for the electric potential, we 
choose to recover the well-posedness of this problem by introducing a regularization of 
the quasi-neutrality constraint. This is the subject of the second part of section 4. In 
section 5, we present a fully-discrete finite volume scheme based on the AP scheme for the 
Euler-Lorentz model coupled with the perturbed quasi-neutrality constraint. Finally, in 
section 6, we present some numerical results which have been obtained with this scheme. 

2 The isothermal two-fluid Euler-Lorentz model. 

2.1 Scaling. 

In this paragraph, we present the scaling of Euler-Lorentz equations which leads to the 
dimensionless following model: 




(2.1b) 



(2.1a) 



n T V x f + q^B], 



e {i,e} , 



(2.1c) 
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where r is the ratio between the ion gyro-period and the characteristic time scale but also 
the square value of the ion Mach number, Tj = 1 and T e are the dimensionleass ion and 
electron temperatures, and e a and q a are defined by 

1 , if a = i, / 1 , if a = i, . * 

q*H _i ifa _„ (2-2) 



e , if a = e, 



where e is the ratio between the unit electron mass and the unit ion mass. Finally, n T , q[, 
qg, r and B correspond to the dimensionless ion and electron density, the dimensionless 
ion momentum, the dimensionless electron momentum, the dimensionless electric poten- 
tial and the external magnetic field and are functions of the position x e IR 3 and of the 
time t > 0. 



Our starting point is the two-fluid Euler-Lorentz model describing a mixture of an ion 
gas and an electron gas. This model writes 

d t n a + V x • q Q = , 

/q a ®q a \ 1 e 
d t (ia + V x • H V x p a = q a (n a E + q a x B) , 

V n a / m a m a 



/( 

^e Q + V x • ( - 

, a e {i,e} 



-q c 



(2.3) 



q« e E • q 



and the ion-electron coupling is insured by the following quasi-neutrality constraint 

ni = n e = n. (2.4) 

In this two-fluid model, rii, q«, e« and (resp. n e , q e , e e and p e ) are respectively the 
density, the momentum, the total energy per mass unit and the pressure for the ion (resp. 
electron) gas. The physical constants rrii, m e , and e stand for the unit ion mass, the unit 
electron mass and the absolute value of unit electron charge. The electric and magnetic 
fields are denoted with E and B and we assume that B is given whereas E is generated 
by the ions and the electrons through the constraint (12. 4|h 

In the present context, we assume that both ion and electron gases are isothermal, i.e. 
we assume that and p e are given by 



Pi = k B Ti rii 



p e = k B T e n e 



(2.5) 



with constant temperatures Ti and T e , and we also assume that the electric field E derives 
from a potential 0, i.e. E = — V x 0- 

Consequently, the model (I2.3p - fl2.4l) is reduced to 

d t n + V x • q Q = , 

k B T a 



d t c[ a + V x ■ + V x n 

\ n J m n 



(2.6) 



k a G {i, e} 



-n V x + q a xB) 
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We introduce characteristic length x, time t, momentum q, ion temperature T, electric 
potential 0, and magnetic field B such that 

x = xx.' , t = tt' , Ti = T, T e = TT' e , (2.7) 



n(xx',tt') = nn'{x',t'), <f)(xx',tt') = 00'(x',t') 5 

q a (xx.',tt') = q^ a (x',t'), B{xx.',tt') = BB'(x',f), 
We consider the natural ratio 



(2.8) 



x n 

q = — , (2.9) 

and we also assume that the electric and magnetic forces are of the same order, which 
means in terms of characteristic scales that 

qB=^-. 2.10) 

x 

We define the characteristic sound speed c for the ions, the characteristic Mach number 
for the ions M and the characteristic cyclotron frequency ZJ for the ions by 



, kuT q e B 

c=\—, M = =, ZJ= . (2.11) 

irti nc rrii 

Considering a low Mach number regime induces 

M=v^, (2.12) 
with r > small and assuming that the applied magnetic field is strong allows us to take 

tco = -. (2.13) 
r 

Finally, we denote the ratio m e /mj with e and we assume that it is a dimensionless fixed 
constant. Then, removing the primed notations and adding r in exponent, we finally 
obtain the rescaled isothermal two-fluid Euler-Lorentz model writing (12. ip . 

2.2 The limit model 

If r converges to in (12.11) . we formally get the model 

<V + V x -q° =0, (2.14a) 
T a V x n° = q a [-n° V x 0° + q° x B] , (2. 14b) 

a e {i,e} , (2.14c) 
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in which the parallel part of q° and q° are implicit. Indeed, if we separate the parallel 
and perpendicular parts of (12.14bl) for any a, we get 



'd t n° + V x -q° =0, 
T a b • V x n° = -q Q n° b ■ V x 0° . 



(q°; 



1 



IB 



■b x (q a T a V x n° + n° V 



1 -°^°), 



, a G {i, e] 



where qi = b x (q x b), b 



(2.15a) 
(2.15b) 

(2.15c) 

(2.15d) 



B 



and II B II 2 = Bl + A 2 + B 2 We observe in this 

reformulated limit model that (q°)j_ and (q°)± can be algebraically computed from n° 
and 0°, but we do not get any explicit constraint for (q°)|| and (q°)||- One way to answer 
to this difficulty is to couple the limit model (" 12 . 1 5[) with the following equations: 



dt((q°)||)-(^(b®b))q° + (b®b)V> 
1 



n° 



a G {i, e} 



+ lim 

t^o Le„r 



■ (b ® b) (T a V x n T + q Q n r V x r ) 



(2.16) 



These equations are not more than the parallel part of (12.1bj) for any a G {i,e}, with 
t — > 0. However, such a coupling is permitted if we are insured that 

b • (T a V x n T + q a n T V X T ) = 0(e a r) , V«e{i,e}. (2.17) 

These hypotheses are coherent with the parallel part of (12.1bj) if we are insured that 



d t ql + V x • (^^) = 0{1) , V« G {i, e} . 



Furthermore, if we consider the limit r — > of (I2.17p . we obtain 

b • (T a V x n° + q a n° V x 0°) = , Vae{i,e} 
which are exactly the equations (12.15bj) . 



(2.18) 



(2.19) 



3 Reformulation of the r-dependent model and of the 
limit model. 



In order to validate the coupling between fl 2 . 1 5 [) and f !2. 16j) . we have to insure that 

b ■ (T a V x n T + q a n T V x r ) = 0(e a r) , VaG{z,e}. (3.1) 
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These properties are validated since the Euler-Lorentz equations (12.11) are equivalent to 



«9 t V V x • ((b <g> b) (T a V x n T + q a n T V X T )) 

= V x • ((ft(b ® b)) «£ - (b ® b) V x • (^^) - 

5t((q r Jii)-(^(b®b)) q ; + (b®b)v x q; 



(3.2a) 



(3.2b) 



+ (b <g> b) (T a V x ri r + q a n T V x r ) = 



<£)± = ^ibx(q a T a V x n T + n T V X T ) 



Ml + v> 



; a 6 {i, e} 



(3.2c) 



(3.2d) 



and more precisely thanks to the diffusion equations (13. 2 aft for both a = i and a = e. 
These equations are obtained for (12.11) by a differentiation in time and position procedure 
and projections in the direction of b and perpendicularly to b. More details about these 
computations can be found in Appendix lAl 
When t — > 0, we obtain 



f V x -((b®b) (T Q V x n° + q a n°V x )) =0, 

<M(q°)||)-(<9*(b®b))q° + (b®b)V 
1 



q° ®q° 



TV' 



+ lim 

T-+0 



t n T 



(b <g> b) (T a V x n T + q Q n T V X T ) 



0, 



(3.3) 



(q 



IB 



-b x (q a T a V x n° + n° V x 0°) 



which is equivalent to (12.1 5ft - ( 127L6]) . 



4 Semi-discrete AP schemes. 

In this section, we propose a semi- discretization of (12. ip in time which is also consistent 
with the reformulated model ( 13. 2p : proceeding in such a way insures us that the approx- 
imation which will be computed ought to this numerical method will be also consistent 
with (13. 3p and, equivalently, with (I2.14p when r converges to 0. 

More precisely, the semi-discretization we describe in the the next lines is based on 
semi-implicit mass fluxes and fully implicit pressure and Lorentz forces. This strategy is 
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motivated by the fact that we want to preserve the balance between the pressure gradient 
and the Lorentz term, i.e. to insure that, at every time step t m , 



T a V x n T ' m + q a [n T ' m V x r '' 



q5 m x 



B m ]=C(e a r), V«G{i,e} 



(4.1) 



where the notation 9 m stands for an approximation of the function 6 = 0(x, t) at the 
time step t = t m . This methodology differs from the AP semi-discretizations which were 
described in [5] and [16] : indeed, in these papers, the authors considered a fully explicit 
mass flux, a fully implicit Lorentz term and a semi-implicit pressure gradient, and this 
leads to a non-conservative discretization of the velocity equation. 

Firstly, we describe the semi-discretization and we reformulate the semi-discrete model 
which is obtained by following the same approach as in Section [3j Then we discuss the 
difficulties which are brought by this reformulation and we introduce a regularization of 
the mass conservation equations ( 12. lap which allows us to bypass these difficulties. 



4.1 Time semi-discretization. 
4.1.1 Asymptotic-Preserving property 

The considered time semi- discretization is the following: 



At 



+ V x • ((b m+1 <g> b m+1 ) q T a m+1 ) 

+ V x ■ ((I - b m+1 ® b m+1 ; 



<fi m ) = o 



q, 



T,m+1 



At 



n T ' m /n T ' m 6?) n r > m \ T 

^- + V x • r° q ° ) + — V x n- m+1 



n 



r,m+l y 



T,m+1 i 



qr.m+l x ,j 



k a e {i,e} . 



(4.2a) 



(4.2b) 



(4.2c) 



As it has been announced above, we chose to implicit the whole pressure and Lorentz 
terms in order to have ( 14. ip at every time step. The choice of the implicitation of the 
parallel part of the mass fluxes is motivated by the fact that we have to reformulate the 
model (14. 2 p by injecting the parallel part of ( 14.2bp with a — i (resp. a = e) in (I4.2ap with 
a = i (resp. a = e). Such a procedure leads to the separate computation of (q[' m+1 )jf +1 , 



,q 



T,m+l\m+l 



^•q-r.m+iw+i anc j (qp m+1 )™ +1 by applying projection operators in the b m+1 
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and orthogonal to b m+ directions. This leads to 



( T,m+l\m+l 
I Ha )± 



q a e a r 
At \\B m+1 :i 

1 

- m+l 



■b m+1 x (q^ m+1 )™ +1 



gm+l I 



b m+i x [q Q T a V x n T - m+1 + n T ' m+1 V x </> r ' m+1 ] 



(4.3a) 



+ 



Igm+l I 



X 



T.m /OS rtTim 



At x V n T - m 



f _T,m+l\m+l 
l^a J II 



(b m+1 ® b m+1 ) q^' m - At (b m+1 g> b 



m+l\ 



v x - 



q r,m q r, m 



At 

k a G {i,e} , 



;b m+1 <8) b m+1 ) [T a V x n T ' m+1 + q a n T ' m+1 V y 



TV 



ir,m+l] 



(4.3b) 



(4.3c) 



on one hand, and to a couple of anisotropic diffusion equations for n T,m+1 and r > m+1 with 
an anisotropy carried by b m+1 (see the computations of Appendix |B] with d = C e = 0) 
on the other hand. These diffusion equations are of the form 



- V x - ((b m+1 ®b m+1 )V> 



n 



r,m+l 



) + \Tn T ' m+1 = tK 



r,m+l 



- V x • (n r ' m+1 (b m+1 <g> b m+1 ) V x T ' m+1 ) = r S T ' m+1 , 
where A only depends on e, At and T e , and where 



K[/\t,l e ,e,n ,q { , q e , b ) 



S(At,T e 



e, n r ' m+1 , n T ' m , q[ ,m , q^' m , b m+1 ) . 



(4.4) 
(4.5) 

(4.6) 
(4.7) 



Remark that the algebraic equations (14.3aj) can be solved for any value of r. Then the 
scheme (14. 2 p is Asymptotic-Preserving if and only if 

b m+1 • (T a V x n T ' m+1 + q a n T ' m+1 V x T ' m+1 ) = 0(e a r) , V a G {i, e} . (4.8) 

These hypotheses are validated if we take into account the following boundary conditions 



(b m+1 ■ V x n r ' m+1 ) (b m+1 • v) = , on dfl, 
(b m+1 ■ V x r ' m+1 ) (b m+1 • v) = , on dQ, 



(4.9a) 
(4.9b) 



alongwith the diffusion equations (I4.4p and ( 14. 5p . Indeed, the solutions of the problems 
f O|) -f l4T9ij) and f^5|) -( l49T5|l satisfy 



jm+1 



• V x ri 



T,m+1 



C(r) 



,m+l 



• V, 



LT,m+l 



0(r) 



(4.10) 



which is ( 14. 8 p up to some linear combinations. 
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4.1.2 Anisotropic diffusion problems 

Now we are insured that the semi-discrete scheme (14.21) is Asymptotic-Preserving, we fo- 
cus on the anisotropic diffusion problems which are satisfied by n T ' m+1 and T,m+1 . 

We remark the diffusion equation (14.4R coupled with the Neumann boundary condition 
(I4.9aj) is well posed for any r > but becomes ill-posed when r = 0. To be more precise, 
the limit of (Q]l -( Oa|l writes 

J -V x • ((b m+1 g> b m+1 ) V x n°' m+1 ) = , on fi, 

\ (b m+1 • V x n°' m+1 ) (b m+1 • v) = , on 90, 1 j 

and a solution n°' m+1 of (14. lip is defined up to a function c : O — > R such that b m+1 - V x c = 
0. Since we want to compute the particular solution n°' m+1 of (14. lip which is exactly the 
limit of (n T ' m+1 ) T>0 when r — > 0, we follow the same approach as in [5] and we use the 
following theorem: 

Theorem 1 (Brull, Degond, Deluzet [5]). Let us consider the subspace K C L 2 (Q) defined 
by 

K= {ueL 2 (n) : b m+1 -V x w = 0}, (4.12) 
and the functional space Wo defined by 

Wo = {u E L 2 (n) : V x • (b m+1 u) E L 2 (n) , (b m+1 • v) u = on dtl } , (4.13) 

provided with the norm ||w||w = ||V X • (b m+1 w)|| i2 (Q\- Then we have the following prop- 
erties: 

1. K is a closed subset in L 2 (Q), 

2. Wo is a Hilbert space and V x ■ (b m+1 Wo) is a closed subset of L 2 (Q), 

3. L 2 (Q) = K®K ± with K ± = V X - (b m+1 W ). 

Having these results in hand and assuming that n T,m+1 is in L 2 (Q), we write 

n r,m+l = n T,m+l + q r,m+l ? ^ ^ 

with 7r T ' m+1 G K and q T ' m+1 G K ± . Since the solution n T ' m+1 of fl4T4D- fl4T9al) is unique 
when r > 0, the functions 7r r,m+1 and q T ' m+1 are also unique as the projection of n T,m+1 
on K and respectively. Then the diffusion problem (14 . 4 j) - f )4T9~a|) writes 

-V x • ((b m+1 ®b m+1 ) V x g r ' m+1 ) 

on O 

+ \ T (7r T > m+1 + q T > m+1 ) = T R T > m+1 , ' (4.15) 

(b m+1 • V x g r ' m+1 ) (b m+1 • v) = , on dVl. 
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If we consider the variational formulation of ( I4.15P over K, we find that (A 7r r,m+1 — 
R T ' m+1 ) e K^, i.e. there exists h T > m+1 such that 

A 7r r ' m+1 - R T ' rn+1 = V x • (b m+1 h T ' rn+1 ) , on tt, 
(b m+1 • v) h T > m+1 = , on dtt. { ) 

By applying the operator b m+1 • V x on this equation, we find that h T ' m+1 is the unique 
solution of 

_ b m+i . v x (V x • (b m+1 h T > m+1 )) = jb m+1 ■ V x R T ' m+1 , on SI, ^ ^ 

(b m+1 • v) h T ' m+1 = , ' on dtt, 

which is a well-posed problem for any r > 0. 

Since q T ' m+1 e K , we claim that there exists / T ' m+1 e L 2 (Q) such that 



q T,m+l = y x . ( b m+l r ,m+l) ^ Qn ^ 

(h m+1 ■ u) l T > m+1 = , on on. 



(4.18) 



(4.20) 



If we consider now the variational formulation of (I4.15P over K ± , we find that / T ' m+1 is the 
solution of a fourth-order problem which can be written as two successive second-order 
problems of the form 

-b m+1 ■ V X (V X • (b m+1 L T ' m+1 )) + t A L r ' m+1 

= _ rb m + 1 .V x ^ m+1 , ° nQ ' (4.19) 
(b m+1 • v) L T ' m+1 = , on dtt, 

-h m+1 ■ V X (V X • (b m+1 r> m+1 )) = L T > m+1 , on tt, 

(b m+1 • v) r> m+1 = o , on an. 

Remark that these problems remain well-posed for any value of r > 0. Then, instead of 
solving the problem f l4.4p - fl4.9ap which becomes ill-posed when r = 0, we solve the prob- 
lems fETTTp . fT4TT9|) and flPU]) for computing ft, r ' m + x and / r - m+1 , then we compute 7r T ' m+1 
and g T ' m+1 by using ( I4.16P and (14.181) respectively, and we finally get n T,m+1 as the sum 
of 7r T ' m+1 and q T > m + l . 

Concerning the problem (I4.5p - fl4.9bp for the electric potential T,m+1 , we remark that 
it remains ill-posed for any value of r > 0. Indeed, assuming that this diffusion problem 
admits at least one solution <f) T,m+1 , we can prove that this solution is not unique: for this 
purpose, we consider a function c : tt — > M. satisfying 



_m+l 



•V x c = 0, ontt. (4.21) 



Then, it is straightforward that T ' m+1 + c is also a solution of the problem f l4.5p -f l4T9bl) . 
Since this proof works for the problem (I4.5p - fl4.9bp and for any value of r > 0, the 
decomposition of its solution does not provide unique projections on K and K x just as it 
is done for n r,m+1 . Then, we have to find a way to restore the uniqueness of the solution of 
the diffusion problem for the electric potential. This is what we do in the next paragraph. 
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4.2 Regularized two-fluid Euler-Lorentz model. 

As it is explained in the previous paragraph, a classical Asymptotic-Preserving scheme 
based on the model ( 12.11) . i.e. based on making implicit the parallel mass fluxes and 
Lorentz and pressure terms, leads to a non-unique solution problem for computing the 
electric potential at time step t m+1 . In order to bypass the difficulty and to restore the 
well-posedness of the problem in T ' m+1 , we choose to include a small regularization in 
the mass conservation equations ( 12. lap . That is why we introduce the terms Cj d t <f) and 
C e dt<p in such a way that this new model writes 



( d t n T + C a d t 4> T + V x -q T a = 0, 



k a e {i,e} . 



+ T a V*n T 

= q a [-n T V x T + q^ x B] , 



(4.22) 



Here Ci,C e > are two fixed small parameters which will be chosen later. Then we 
consider the same semi-discretization method as previously, i.e. 



n 



T,m+1 



n 



LT,m+l 



At 



At 



+ V x • ((b m+1 ® b m+1 ) q T a m+1 ) 

+ V x • ((I - b m+1 <g> b m+1 ) q£ m ) = , 



(4.23a) 



At 



+ V x - 

la 



H V x n 



T,m+1 



[ 



n 1 '"' / e a t 
n T ' rn+1 V x </> T ' m+1 



+ cQ m+l x B m+1 ] 



(4.23b) 
(4.23c) 



, a G {«, e} . 

By splitting parallel and perpendicular parts of ( 14. 2 3b j) for a = z and a = e, we get 
that (qa' m+1 )™ +1 and (qa' m+1 )jj 1+1 satisfy (I4.3p . Following the same procedure as in the 
previous paragraph (see Appendix [B]) , we inject (14.3bl) with a = i (resp. a = e) in ( I4.23al) 
with a = i (resp. a = e). Under the hypotheses 



d + eC e = and Q - —C e = C . 



(4.24) 



with C > being given, we find that n r,m+1 and r > m+1 satisfy the following diffusion 
equations: 

- V x ■ ((b m+1 ® b m+1 ) V x n r ' m+1 ) +rAi n T ' m+1 = r _R T > m+1 , (4.25) 

- V x ■ (n T ' rn+1 (b m+1 g> b m+1 ) V x T ' m+1 ) + r A 2 r ' m+1 = r S T ' m+1 , (4.26) 
where Ai, A 2 only depend on e, At, T e and C, and where 



T,m+1 



R{At,T e 



e, n 



r,m „ T > m i -, r > m '^ m +l" , 



(4.27) 
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S r,m+1 = ^ J At> £j n r, m+ l ; n r, m? ^ q r,m^ ^ ^+1) _ ( 42g ) 

Remark that the constraints (I4.24p are equivalent to 

T e C T e C 

Ci = ^—, C e = , — . 4.29 

1 + T e ' e(l+T e ) v 7 

As a consequence, it is necessary to take C > small enough to insure that Cj and C e 
are close to 0. For that, we can take C = 0(e) provided that the ratio e = m e /mj is small. 

We couple (14.251) with the boundary condition given in (I4.9al) . This diffusion problem 
is well-posed for any r > and becomes ill-posed if r = because of a lack of uniqueness 
of the solution (see page[TU]). As a consequence, we can apply Theorem [T] and write n T ' m+1 
under the following form: 

n T,m+l = n r, m +l + q r,m+l ^ (4 3^ 

with vr r ' m+1 G K, q T > m+l G K x defined by 



7T 



T,m+1 



1 



Ai 1 x v 7J ' (4.31) 

where h T,m+1 and / T ' m+1 are the solutions of 

-b m+1 • V X (V X • (b m+1 h T ' m+1 )) = yb m+1 • V x /T' m+1 , on n, ^ 32 ^ 
(b m+1 • v) h T ' m+1 = , on dQ, 



and 



-b m+1 • V X (V X • (b m+1 L T ' m+1 )) + t Ai L T ' m+1 

= -rb m+1 .V x iT'' w+1 , ° n ^' (4.33) 
(b m+1 • 1/) L T ' m+1 = , on <9ft, 

-b m+1 • V X (V X • (b m+1 Z T ' m+1 )) = L T > m+1 , on fi, 
fb m+1 • v) l T > m+1 = , on 0ft. 



(4.34) 



Now, we couple (14. 26ft with the boundary condition given in fl4~9bl) . This diffusion 
problem has the same properties of (I4.25j) - (l4.9aj) which has been discussed above. Then, 
we can write T ' m+1 under the following form: 

0T,m+l = ~r,m+l + ~r,m+l ^ ^35) 

with 7f r ' m+1 G K and <f ' m+1 defined by 



~T,m+l 



1 



7T 



j^m+l + a 2 V x ■ (b m+1 /T' m+1 )] , 



g r ' m+1 = V x - (b m+1 f r ' m+1 ), 
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where h T,m+l and l T > m+1 are the solutions of 

-b m+1 • V X (V X ■ (b m+1 h T ' m+1 )) = yb m+1 ■ V x 5 T ' m+1 , on ft, ^ 
(b m+1 ■ v) ~h r ' m+1 = , " on <9ft, 



and 



-b m+1 • V X (V X - (b m+1 L T ' m+1 )) +T\ 2 L T ' m+1 

on 

= -rb m+1 -V x 5' r ' m+1 , ' (4.38) 

(b m+1 ■ u) L T > m+1 = , on <9ft, 



-b m+1 ■ V X (V X ■ (b m+1 l T > m+1 )) = L T > m+1 , on ft, 
(b m+1 • z/) l T ' m+l = , on <9ft. 



(4.39) 



5 Fully-discrete scheme 



In this section, we present the fully-discrete version of the Asymptotic-Preserving method 
for (I4.22p we have presented in the previous paragraph. Before going further, we introduce 
some notations which will be used throughout this section. 

First, we consider a uniform mesh (xj, yj, Zk) = (i Ax, j Ay, k Az) on ft and we define the 
following subsets of Z 3 : 



I = {(h3,k) G Z 3 : (x u yj,z k ) E ft} , 

I = {(i + a,j + P,k + -r) : (i,j,k)el, (a,/3, 7 ) G {-1,0,1} 3 }. 
Then, we define the meshed domain ft^ by 



(5.1) 



^h= [J [Xi-i/2, X i+1 / 2 ] X [yj-i/2,yj+l/2] X [Zk-1/2, Z k+ l/ 2 ) ■ (5.2) 

(i,j,k) £ / 



We also define the subsets and 7* of Z 3 as follows: 

7* = {(«, j, fe) G Z 3 : (^+1/2,^+1/2,^+1/2) G ft/J , 
7* = {(«, j, k) G J* : (x i+ i/2,yj+i/2,Zk+i/2) <£ dQ h ] . 



(5.3) 



Finally, we assume that, for any K = (K x , K y , K z ) G Z 3 , the notation "\k" stands for 
an approximation at the point (xK x ,UK y , zk z ) and that the notation "\k* v stands for an 
approximation at the point (xk x +i/2, UK y +i/2, Zk z +i/2)- From now, we also denote the 
point {x Kx ,y Ky ,z K J by a cell center and the point (x Kx+1/2 ,yK y +i/2, z Kz+1/2 ) by a node. 

The next lines are structured as follows: firstly, we present the finite volume scheme 
based on the semi-discretization (I4.23p . Then, we reformulate the obtained fully-discrete 
scheme by following the same approach as in section 4 and we suggest a numerical method 
for solving the fully-discrete diffusion equations for n T,m+1 and r,m+1 . 
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5.1 Finite volume scheme 

First, we introduce some notations for the explicit and implicit fluxes for the hydrody- 
namic part of 1 14.231) : 



eexp,T,m 
a,a 



I e • ((I — b m+1 <g> b m+1 ) c& m ) \ 

/ 



(5.4) 



n' 



/ e a ■ ((b m+1 ® b m+1 ) q^ m+1 ) \ 



r,m+l 



V 



(5.5) 



/ 



where a G {«, e}, a G {x, y, 2}, e a is a vector of the canonical basis of M. 3 , and where I is the 
3x3 identity matrix. Then, we consider different notations for divergence and gradient 
operator depending on whether they are applied on some component of the implicit or 
the explicit fluxes. More precisely, we define the operators Vf v -, VV and by linking 
them to the fluxes by the following relations: 



V\ F ■ ((I — b m+1 ® b m+1 ) cQ m ) 

^ V n T ' m 
V h ■ ((b m+1 ® b m+1 ) q^ m+1 ) 



r,m+l 



Ea rimp,T,m+l 

a £ {a;,2y,2} 



(5.6) 



(5.7) 



where a G {«,e}. As a consequence, the fully-discrete model obtained from (I4.23[) is 
written as follows: 



n T ' m+1 \ 



n T ' m \ 



At 



+ (y h ■ ((b m+1 ® b m+1 ) q; ,m+1 )) k 
+ (v£ y ■ ((I — b m+1 <g> b m+1 ) q ;- m ) 



(5.8a) 



-a 



tTjm+l d) T ' m 

\k i_ \k 

At 



At 



/ a T ' m (X) a T ' m 



1 



k a G {i,e} . 
In finite volume terms, we have 

\ U a'-a,a )\k ~ 



[T a V h n T ' m+1 + q a (n T ' m+1 V h ^ m+1 - c£ m+1 x B m+1 )] 



Aa 



I -pr,m 

a,a 



\K+e a /i 



-T. 



a, a 



K-e a /2' 



(5.8b) 
(5.8c) 

(5.9) 
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with 



7 T 



I cexp,T,m I cexp,T,m \ 



-© r ' m , (W T ' m , - w r - m 



lif+ea ' 
Ik+oo 



(5.10) 



where a G {z, e} and a G {x, y, z}. In these formulae, W^ m | K is 



n r,m. 



and D^,'™ is the numerical viscosity matrix linked with the flux f^P' 1 "'" 1 for any a G {i, e}. 
Since the main goal of the present paper is to validate the time semi-discretization pre- 
sented in the paragraph 4.2, we choose to compute the viscosity matrices with Rusanov's 
method (see [3D] and [3D])- For this purpose, we denote the eigenvalues of the jacobian 
matrices Jac Wa (f«T T ' m ) by A££ a (k = 1,2,3,4). Then, D£™ is defined by 



a,a 



K + e a /2 



I max max (|AT;? al | , |A^T al |). (5-12) 

fc=l, 2,3,4 V a ' k ' a \K+ ea 1 1 a : fc . a | K 7 v ' 



5.2 Reformulation of the fully-discrete scheme 

Following the same approach as in sections 3 and 4, we reformulate the discretized model 
( 15. 8j) by computing separately the perpendicular part and the parallel part of q[' m+1 | K 
and qe' m+1 | x and by solving two diffusion equations to find n T ' m+1 \ K and 4> T,m+1 \ K - More 
precisely, under the hypotheses (14.24)) for Cj and C e , we solve 

- (v h -((V h n T > m+ X +1 )) +\ 1 m T ' m+1 lK = TR T ' m+1 lK , (5.13) 

^ ' \k 

for finding n T ' m+1 \ K and 

- (V h • K< m+1 (V h T ' m+ \ +1 )) + A 2 r<^ m+1 k = rS T ' m+1 llc , (5.14) 

^ ' I K 

for finding r ' m+1 | K - These discrete diffusion equations are obtained by injecting the 
parallel part of (15.8b)) with a = i (resp. a = e) according to h m+1 \ K into (15.8a)) with a — i 
(resp. a = e), then performing some linear combinations of the obtained equations (see 
Appendix O) . 

Having n T,m+1 | A , and <f) T,m+1 i K in hand, we can compute separately the parallel part 
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and the perpendicular part of q[' m+1 | x and ql ,m+1 \ K by using the following formulae: 

ll<la )\\ )\ K 

((r„v^"" +1 + ( , t ,^™ +1 v^"' +1 ); +1 ) 

\ II / \ K 



(5.15a) 



At 



((q[' m+1 )T +1 ) 



q " ' " 7 b m+1 k x ((q[' m+1 )r +1 ) 



lx At||B m+1 | A . 



1 



. ||B m+1 
q Q e„r 



b m+l x ( n r,m+l Vh ^' m+1 + q a T a V h U T ' m+1 ) 



+ 



B 



m+l 



O r ' m , / n T ' m « n T ' m 



n 1 



k a G {i, e} 



(5.15b) 



(5.15c) 



5.3 Three-point scheme 



In this paragraph, we focus on the resolution of ( 15. 13)) and (15 . 14)) provided with a dis- 
cretization of the Neumann-like boundary conditions (14. 9p . For solving these diffusion 
equations, we follow the approach of Degond & Tang in [19]: we choose a three-point 
scheme by replacing the equations ( 15 . 13)) and (15. 14() by 



(d™+V' m+1 )| 



T,m+1 



tK 



r,m+l 



Vif el, 
VKeZ\l* 



(5.16) 



and 



■(^ +1 «' m+1 ^r + V T ' m+1 )) , +r A 2 



r,m+l 



(srv^ 1 )^. = o 



QT,m+l ' 

T D |k 5 _ 

\/Keh\h 



(5.17) 



respectively. In these equations, ^*' m+1 | K stands for the average of n T,m+l on the node 
(xk x +i/2, VKy+i/2, z Kz+1/2 ) defined by 



n 



r,m+l 



m+l 



(5.1* 



a,/3, 7 e{0,l} 
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and the operators d™ +1 and d™^ 1 correspond to some approximations of b m+1 • V x and 
V x • (b m+1 -) respectively. These operators are defined by 



(fle? 1 p)k= E 

/3,7 6{0,1} 

+ E 

a,7£{0,l} 

+ E 

a,/3e{0,l} 



/ E 

/3,7e{0,l} 

E 

a, 7 G{0,l} 

E 

\ «,/3 6{0,l} 







P\K + pe y + 1 e z 




AAx 






+ By + ~fe Z 


P\ K + ae^ + 7e 2 




4Ay 




^ \K + ae x 


+ pe y + e z 


'P\K+aLex+(3ey 



\ 



AAz 



(5.19) 



A'* — e^ — pe y — ~fe 







AAx 


{tip 




- (b m+1 v)\ 






AAy 




- l P)\ K , 


- (b m+1 v)\ 

-aex—pey V Z ± / | Km —ae x — ,8e y — e z 



(5.20) 



AAz 



Then, replacing b m+1 ■ V x and V x • (b m+1 -) by <9™ +1 and d™^ 1 in the decomposition 
procedure (I4.30l) - fl4.39p . we compute n T ' m+1 \ K and (fi T ' m+1 \ K that satisfy 

(T a ^ + V- m+1 + q a <' m+1 ^ + V r ' m+1 )k, =0(e a r), \/ae{t,e}. (5.21) 

We remark that we have the good Asymptotic-Preserving property on the nodes. However, 
we need it on the cell centers, i.e. 



(b m+1 ■ (T a V h n^ m+1 + q a n T ' m+1 V h ^ m+1 )) ]r = 0(e a r) 



(5.22) 



for any a G {i, e}. To reach such a result, we introduce the discrete gradient * defined 
by 



E 

/3,7e{o,i} 

E 

a,7£{0,l} 

E 

\ «,/3G{0,l} 



P|x+e a; + / 3ej / + 7e z 


P\K+0e y +ye z 


4Ax 




PlK+OLex+ey + ryez 


V\ 


AAy 




P\K + ae x +0 ey + ez ~ 


P\K+ae x +fley 


AAz 



(5.23) 



and we use it to couple the three-point scheme we have presented with the formulae for 
((q? W+1 )T +1 ) k) ((qe' m+1 )r +1 ) k; ((qr +1 )[T +1 )| A . and H' m+1 )\\ +1 )\ K : more precisely, 
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we use the operator Vh,* for computing the following terms in (15.1 5ft : 



Ik- 



T,m+1 



a,^,7S{0,l} 



r,m+l 



7 iT,m+l\m+l\ 
h,*<P J|| J 



X* — ae s - /3e y — 7e 



( ||B"+'ii b "' + ' x i - q " nT ' m+1 V ^ r,m+1 + r«v t ^ m+1 ; 



- E 



a,/3,7S{0,l} 



|B m+1 



_ h m+l x ^r,m+l y^^r.m+l 



(5.24a) 



B X —pey—"fe z 



(5.24b) 



(5.24c) 



As a consequence, we obtain the properties (I5.13P and (15.141) on cell centers and we 
can compute the parallel part and the perpendicular part of q[' m+1 and q^> m+1 by using 
separately the formulae (I5.15p . 



6 Numerical results 

In this last section, we present some 2D numerical results which have been obtained 
with the AP scheme we have presented in sections 4.2 and 5 for the perturbed two-fluid 
euler-Lorentz model (14.221) . 

6.1 Validation of the three-point scheme for the diffusion prob- 
lems 

Since the AP scheme relies on the truthfulness of the properties (14.81) . we first present 
some numerical results from the three-point scheme used for the resolution the diffusion 
problems for n T ' m+1 and r ' m+1 (see paragraph 5.3). More precisely, the main goal of the 
first test sequence is to solve the diffusion problems (l4.25p -( l4T9al) and (I4.26p - (l4.9bl) for any 
value of t > and to insure that 

WK, n T ' m+1 lK ^n°' m+1 lK , <p T ' m+1 \ K -> <P°' m+1 \ K , (6.1) 

as t — > 0. To perform this validation, we apply our method to the following diffusion 
problem: 

f -V x - (^(b®b)V x p r ) + rA^ = rf , on ft, 

\ (h t (b ® b) • v x p r ) • v = o , on an, 1 ' 
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where A > 0, f T : Q -> E, if T : Q -> and b : ft -> M 3 are given. 
Let us consider a function sequence (p T ) r >o denned by 

P T =Po + rp T 1 , 

with p and p\ satisfying 



and 



b ■ V x po = , on O, 

(h t (b ® b) Vxpj) • v = o , on an. 



We assume from now that 

f = Ap T -V x - (# r (b®b)V x K), 
which implies that p T is the analytic solution of the problem (16. 1\ . 



(6.3) 

(6.4) 
(6.5) 

(6.6) 





Figure 1: L , L and L°° norms of the error between p T and its approximation p T app as 
functions of h: case with r = 10~ 2 (left) and r = 1CT 9 (right). 

In Figure HJ we plot the evolution of the error between p T and its approximation 
(denoted with p T app ) as a function of the space step h. In these results which are presented 
in decimal logarithmic scale, we have chosen r = 10~ 2 and r = 10 -9 . Q is set to [1,2] x 
[1,2] C M 2 , p and p\ to 



Po{x,y) =2, Pi(x,y) = ({x - 1)(2 - x)(y - 1)(2 - y))* 
and A, if r , and b are chosen as 

A = l, 

H T (x,y) = 1 + sin 2 (x) sm 2 (y) , 
b = (sin#, — cos 9) , with 8(x, y) = arctan(y/x). 



(6.7) 

(6.8) 
(6.9) 
(6.10) 
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Figure 2: L 1 , L 2 and L°° norms of the error between p T and p° = po as functions of r: 
case with a 100 x 100 uniform mesh. 



As we can remark in this figure, the solver for the diffusion problem ( 16. 2j) based on 
the micro-macro decomposition presented in Section 4.1 and on the discrete differential 
operators dh and dh,* is second order accurate in h since we observe that the error \\p T — 
PappllLP ip — 1> 2, oo) is linearly decreasing in log 10 scale when h — > 0, with a slope which is 
equal to 2. This is due to the fact that, according to the definitions (I5.19P and ( I5.20p . the 
operators dh and dh,* are themselves second order accurate. Furthermore, we have this 
property for r = 10~ 2 and r = 10~ 9 , so we can conclude that the second order accurate 
of the solver is not penalized by the smallness of r. 

Together with this convergence results in h, we plot in Figure |2] the error between p T app 
provided by the solver and p° = po as a function of r with a 100 x 100 uniform mesh, and 
we take the same values of po, p\, A, H T , and b as above. By definition of the analytic 
solution of p T and po (see ( 16. 3p ). we except this error to be of the same order of r. This is 
confirmed by FigureEJ indeed, the error \\pl pp — Po\\lp (p = 1, 2, oo) is linearly decreasing 
in log 10 scale when r converges to with a slope which is equal to 1. 

From these two results, we can claim that 

limlirn^ =Po, (6.11) 
h — v t y ^ 

which is exactly to say that the numerical solver for the diffusion problem (16. 2p is 
Asymptotic-Preserving when r — > 0. Then, we can use it for solving the diffusion problems 
(OoP - dPa]) and P^ - lpHl for n T > m+1 and r - m+1 . 



6.2 Numerical results for the two- fluid Euler-Lorentz model near 
the drift-fluid limit 

In order to validate the AP scheme we have developed for the perturbed Euler-Lorentz 
(14.221) . we compare the results which are computed by the AP scheme to those which can 
be produced with a fully explicit finite volume method. From now, we denote with classical 
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method a finite volume scheme which is based on the following time semi-discretization: 
At +C « ' At +V X .^ = 0, 



At " x V n^™ / e a r x (6.12) 



q a \ n T ' m V x T ' m + q£ m+1 x B m+1 l 



a E {i, e] 



We easily remark that the stability condition of such a method strongly depends on r: 
as in a low Mach number numerical experiment, the smaller r is, the smaller the time 
step At must be in order to insure that a method based on (16.121) is stable. As a ex- 
ample, if we solve the hydrodynamic part of (16.121) with Rusanov' scheme, we must have 
At = Oihr 1 ! 2 ') at least, where h = min(Ax, Ay, Az). 

In the next lines, we distinguish two opposite situations: 

• The resolved case: The time step is small enough in order to insure that both 
classical and AP methods capture the fast time variations within the solution, 

• The under-resolved case: The time step does not allow the capture of fast time 
variations but insures at least the stability of the AP scheme. 

The test case we present here is based on the perturbation of the following stationary 
case: 

• The magnetic field is uniform and reads B = (sin a, — cos a, 0) with a6l fixed, 

• n T '°(x,y) = no, 4> T '°(x,y) = </>q, q['° = qg' = B with some constants n and (j>o, 

• The whole system (14.221) does not depend on the variable z. 

Remarking that both classical and AP schemes compute the exact solution provided with 
these initial datas, we choose to introduce a small perturbation at the initial time step. 
More precisely, we choose to replace n T '°(x,y) = n by 

n T '°(x, y) = n + r max (0, 1 - 77 (x - x ) 2 - i] (y - y ) 2 ) , (6.13) 

with 77 > and (x , yo) G 
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qj x (AP scheme) 




qj y (AP scheme) 



Figure 3: Resolved case at time t — 6 x 10 
q[ as functions of (x, y) computed with t 
(right). 




q[ x (classical scheme) 




q T i y (classical scheme) 



: x and y components of the ion momentum 
i AP scheme (left) and the classical scheme 
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q T ex (AP scheme) 



q T ex (classical scheme) 




_ 




ql (AP scheme) 



ql y (classical scheme) 



Figure 4: Resolved case at time t — 6 x 1CP 6 : x and y components of the electron 
momentum ql as functions of (x, y) computed with the AP scheme (left) and the classical 
scheme (right). 



We also assume that the physical domain Q is [1, 2] x [1, 2] and is meshed by a 100 x 100 



uniform mesh. We also precise the initial datas by taking r 



10" 



1 T 



(7=10- 



0. 



rj = 80, (x , yo ) = (§>§)» "o = 1 and ^ 
In Figures [30 we present some results in the resolved situation for both classical and 
AP schemes and At = 5 x 10 -9 is taken as time step. As we can see in these figures, 
the results which are produced by the AP scheme are very close to the classical method's 
ones. 
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qj y (AP scheme) q[ y (classical scheme) 



Figure 5: Under-resolved case at time t = 6 x 10 : x and y components of the ion 
momentum q[ as functions of (x, y) computed with the AP scheme (left) and the classical 
scheme (right). 
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q T e (AP scheme) q T e y (classical scheme) 



Figure 6: Under-resolved case at time t = 6 x 10 : x and y components of the electron 
momentum cf e as functions of (x, y) computed with the AP scheme (left) and the classical 
scheme (right). 



In Figures [5J1S1 we present some simulations which are obtained in the under-resolved 
case, i.e. where the time step is taken much larger than the time step which is required to 
insure the stability of the classical scheme. In the present case, we have chosen At = 10 -6 , 
which is 200 times larger than the time step which has been used for the resolved case 
above. As we can see in these figures, the AP scheme remains stable and produces the 
same results as in the resolved case. However, the classical scheme blows up after a small 
number of time iterations, which is not surprising because the time step we have chosen 
is too large for satisfying the stability condition of this scheme. 

From this numerical experiment, we can conclude that the AP scheme we have de- 
veloped for the perturbed two-fluid Euler-Lorentz model (14.221) allows us to take a time 
step which does not satisfy the stability condition required for capturing the fast time 
variations of the solution. Furthermore, such a time step choice does not penalize the 
quality of the results which are obtained with the AP scheme. 
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6.3 Impact of the perturbation parameter C 

Since the AP scheme we have built for the resolution of the perturbed Euler-Lorentz 
model ( 14.221) has been validated in terms of quality of results, the impact of the value 
of C needs to be investigated. Indeed, we recall that this parameter is linked with the 
constants Ci and C e by the relations 



Ci 



T e C 

1 + T e ' 



e 1 + T e 



(6.14) 



and these constants are introduced to recover the uniqueness of the solution for the diffu- 
sion problem (14.5jl -( f479bl) (see Section 4). Since Ci and C e are introduced in ( 12.11) through 
a perturbation of the mass conservation equations, these constants are assumed to be as 
close to as possible. However, this is equivalent to assume that C is close to 0, so the 

J 1 Q 

diffusion problem (I4.26l) - (l4.9bl) is ill-conditioned since A2 = , r (see Appendix 

At 2 (T e - 1) 

IB"]) and degenerates into the non-unique solution problem ( 14.5jl -( f479bl) when C — > 0. Then, 
it is necessary to investigate the consequences of the choice of C on the stability of the 
AP scheme. 




At = 10 




At = 10 




At = 10" 



Figure 7: qJ tX at time t = 6 x 10~ 6 with C = 10~ 2 and At = 10~ 6 , 10~ 7 , 10 
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In the last test sequence, we run the AP method with the initial datas which have 
been used in the previous paragraph, i.e. 

• We take a 100 x 100 uniform mesh over Vt — [1,2] x [1,2], 

• The magnetic field is uniform and is defined as B = (sin a, — cos a, 0) with a = 

• The initial electric potential is <f) T, ° = <fio with 0o = 0, 

• The initial ion momentum and the initial electron momentum are defined by qp° = 
q!° = B, 

• The initial density is n T, ° defined by 

n T '°(x, y) = n + r max (0, 1 - 77 (x - x ) 2 -r)(y- y ) 2 ) , (6.15) 
with n = 1 constant, r\ = 80 and (xo, yo) = (§, |), 

• We choose T e = 3, e = 1 and r = 10~ 8 . 

In Figure [7J we plot the x-component of the ion momentum q[ which is obtained with 
the AP scheme (I4.23P at time t = 6 X 10~ 6 with C = 10 -2 and a time step At which is 
equal to 10" 6 , 10" 7 and 10" 8 . We can remark that the AP method is stable for At < 10 6 
since all the results are similar. 

Now, we do again this numerical experiment with C = 10~ 3 instead of C = 10 -2 . As 
we can remark in Figure [8] where qj x is plotted at time t = 4 x 10 -6 , the AP method is 
stable with At = 10~ 7 and At = 10~ 8 . However, we observe some important boundary 
effects in the At = 10~ 6 case. 

Finally, we perform this numerical experiment with C = 10~ 4 . In Figure O we plot 
qj x at time t = 2 x 10 -6 with At = 10~ 6 , 10~ 7 , 10~ 8 . No numerical artifacts are present 
in the At = 10~ 8 case whereas we observe some boundary effects if At = 10~ 6 and even 
the blowing up of the method when At = 10~ 7 is considered. 

The results which are provided by the whole test sequence above indicate that we have 
a stability condition for the AP scheme which clearly depends on C and which would be 
of the form 

At = 0(C) . (6.16) 

As a consequence, we can say that the AP scheme we have developed for the perturbed 
two-fluid Euler-Lorentz model (14 . 2 2 [) is Asymptotic-Preserving when C > is fixed and 
when t — > 0. However, it is not Asymptotic-Preserving when C — > and r > is fixed. 

7 Conclusions and perspectives 

In this paper, we studied the isothermal two-fluid Euler-Lorentz system coupled with a 
quasi-neutrality constraint in a low Mach number regime and a strong magnetic field 
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regime. After having presented the model and its limit regime, we proposed a reformula- 
tion of this model which is compatible with the construction of an Asymptotic-Preserving 
scheme. Then we presented the time semi-discrete AP scheme itself and its reformulation 
leading to the resolution of some anisotropic diffusion equations for n T,m+1 and T ' m+1 . 
The equation for T ' m+1 being ill-posed, we restored the uniqueness of the solution of 
this equation by introducing a small perturbation in the mass conservation equations. 
Finally, we performed some numerical tests of this scheme by comparing the results from 
the AP scheme to those from a fully explicit method and we tested the influence of the 
perturbation parameter C on the behaviour of the AP scheme. 

At this point, several work pathes can be investigated. The first one is to go back to 
the time semi-discretization described in section 4.1 by invoking boundary conditions for 
the computation of the electric potential which are different from Neumann conditions. 
The second one is to generalize the present work to some two-fluid Euler-Lorentz models 
involving other pressure laws for pi and p e . 
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At = 10 




At = 10 




At = 10" 



Figure 9: qT x at time t = 2 x 10~ 6 with C = 10" 4 and At = 10~ 6 , 10~ 7 , 10~ 8 . 

A Reformulation of the Euler-Lorentz model 

In this paragraph, we detail the reformulation procedure of the Euler-Lorentz model 

f d t n T + V x -<& = 0, (A.la) 

+T Q V x n T ^ Alb ^ 
= q Q [ - n T V X T + x B] , 

(A.lc) 



5 t q Q + V x ■ 
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leading to the model (I3.2p . First, we separate the parallel and perpendicular parts of 
(lA.lbj) for each value of a: we obtain 



r$n T + V x .(c£)|| + V x -(<£)± = 0, 

M(*)n)-(^( b ® b )) < £ + ( b ®b)v > 



q; ® q; 



(A.2a) 
(A.2b) 



B| 



+ (b ® b) (T a V x n r + q Q n T V x r ) = 

e a r 

-bx (q a T a V x n T + n r V x r ) 



+ W bX 



rr 



. a G {i, e} 



(A.2c) 



(A.2d) 



In order to obtain (13. 2 aft for any a G {i, e}, we compute the divergence in space of ( 1A.2bl) 
on one hand and the derivative in time of ( 1A.2al) on the other hand. We obtain 



V x ■ (<9 t ((qD]|)) - V x • ( (d t {b ® b)) c£ - (b ® b) V, 



q^® q» 



and 



+ V x - ((b®b) (T Q V x n T + q a n T V x T )) =0 



#n T + <9 t V x ■ (q^)|| + 0*V X • (q;)± = 



(A.3) 



(A.4) 



for any a G {i,e}. By doing some linear combinations of these equations, we obtain a 
system of 2 equations for n T and <p T of the form 



<9 t V V x • ((b ® b) (T Q V x n r + q a rf V x r )) 

= V x • ((d t (b ® b)) c£ - (b ® b) V x ■ (^^) " ft((<£)±) 
a G {z,e} , 



(A.5) 



which are exactly the equations (I3.2a|) . Then the Euler-Lorentz model OA.ip is equivalent 
to the combination of ( 1A.2bj) . (1A.2cj) and flA.5j) . 



B Reformulation of the semi-discrete problem 

This paragraph is devoted to the reformulation of the semi-discrete problems (14. 2 j) and 
(I4.22p . Since the model (14. 2 p is not more than (I4.22p with Cj = C e = 0, we present the 
reformulation of the semi-discrete scheme for the perturbed Euler-Lorentz model. This 
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scheme is recalled here: 



( n r,m+l 



n 



LT,m+l 



LT,m 



At 



At 



+ V x • ((b m+1 ® b m+1 ) q T a m+1 ) 

+ V x • ((I - b m+1 <g> b m+1 ) c£ m ) = 



(B.la) 



q, 



T,m+1 



q, 



At 



k a e {z, e} . 



+ V> 



n 1 



+ 



a 



V x n 



T,m+1 



la 



.T,m+1 



[ - n T ' m+1 V x T ' m+1 + cfc m+1 x B m+1 ] , 



(B.lb) 



(B.lc) 



First, we separate the parallel and the perpendicular parts of flB.lbj) according to 
b m+1 . More precisely, we obtain the equations (14.3aj) by performing the vector product 
of b m+1 by flB.lbj) . Concerning the equations (I4.3bj) . we obtain them by multiplying the 
tensor (b m+1 ®b m+1 ) by flB~Ib|) . 

In order to obtain the diffusion equations ( I4.25P and (I4.26p . we put ( I4.3b|) with a = i 
(resp. a = e) in (IB.laj) with a = i (resp. a = e). We obtain a system of two diffusion 
equations satisfied by (n 7 



,T,m+l 0T,m+l x 



' n T,m+l n T,m 



LT,m+l f k'r,m 



q a 55 q a 



At " At 

At V x • f (b m+1 <g> b m+1 ) V x • ( ^-^ 

V V n T ''' 

M V x • ((b m+1 g> b m+1 ) (T Q V x n T ' m+1 

+ q a n r ' m+1 V x r ' m+1 )) = 0, 

k «6{i,e}. 

By doing some linear combinations, these diffusion equations write 
-V x ■ ((b m+1 <g> b m+1 ) V x n T ' m+1 ) 



(B.2) 



+ r 



1 + e , . d + e C e 

n T,m+l _|_ T ' AT,m+l 



At 2 (1 + T e 



At 2 (l + T e ) 



1 + T e 



^V x -(q? +e ^) + _„^ + _^- 



(B.3) 



+ V x ■ (b 



m+l i_,rra+l\ 



<q* ®q* x 



+ eV -"( ^ )J 
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and 



V x • (n T ' m+1 (b m+1 g> b m+1 ) V x T ' m+1 ) 



T P Ci e C. 



T P -e 



At 2 (T e - 1) 



At 2 (T e - 1) 



T,m+1 



T X 



T P -1 



At 



Vx-(qr-^qr) 



At 2 T f 



n 



+ 



TfCi — e C e 
At 2 T. 



T ' m + V x ■ ( (b m+1 ® b m+1 ) [V, 



Mi ®q» v 

n r,m J 



^ v ql' m (8>ql' m 



n' 



(B.4) 



If we consider the constraints (14.241) for Cj and C e , we make the diffusion equations 
uncoupled. Firstly, we compute n T ' m+1 by solving 



V x ■ ((b m+1 ® b m+1 ) V x n r ' m+1 ) + r Ai n T ' m+1 = r iT 



m+l 



(B.5) 



with 



R 



A I 



T,m+1 



1 + e 



At 2 (1 + T e 



1 + T e 



V7 / T.m i T,m\ i ^ T 

" Al Vx ' (q * +eqe )+ A^" 



+ V x - ((b m+1 ®b m+1 ) [V, 



r,m „ r,m 

®q» s 

n T,m J 
^ n r,m 



(B.6) 



then we use it to compute T ' m+1 by solving 



V x • (n T ' m+1 (b m+1 ® b m+1 ) V x r ' m+1 ) + r A 2 r,m+1 = r S T ' m+1 , (B.7) 



33 



with 



A 2 

gT,m-\-l 



T e C 



At 2 (T e - 1) ' 



T e -\ 



— V ■ (a T ' m - — a T ' m ) 
At x [% T^ e > 



-T, 



At 2 T f 



e ( n T ' m+1 - n T ' m ) 



C 



At 2 



+ V x • ( (b m+1 ® b m+1 ) [V, 



(B. 



C Reformulation of the fully-discrete problem 



In this paragraph, we develop the reformulation procedure for the finite volume scheme 
which is detailed in Section 5. This scheme writes 



/ n T,m+l 



n 



At 

- (vr-((i 



^+ (V h - ((b m+1 ® 

b m+l ^ b m+l) q r,m) 



y-m+lN T,m+1 



)), 



K 
T,m+1 



(C.la) 



q r,m+l _ q r,m T 
*" Ik I (V FV ( 



At 



,T,m 



n' 



«£)) 

* K 



E {i, e] 



[T a V h n^ m+1 + q a K' m+1 V^ 



T,m+1 



At 



q r,m+l x B m+1)] 



(C.lb) 



(C.lc) 



As in the previous appendices, we separate the parallel part and the perpendicular 



part of (IC.lbj) according to b m+1 | A ,. By multiplying the tensor (b 



m+l 



b m+ V) b y 



(IC.lbj) . we obtain (I5.15ap for each a. Concerning ((q^' m+1 )™ +1 ) , we compute the vector 
product of b m+1 | K and (IC.lbj) and we obtain (15.15bj) . 

In order to obtain the discrete diffusion equations (15.13j) and (15. 141) . we follow the 
same procedure as in the semi-discrete case (see Appendix [B]): we replace (q[' m+1 )™ +1 
(resp. (ql' m+1 )n l+1 ) by its expression given by (I5.15aj) with a = i (resp. a = e) and we 
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obtain two diffusion equation of the form 



-\\/h ■ {{-La Vhn T + q a n T Vh<P )\\ ) ) 



+ ^K' m+1 + C0 r ' m+1 ) 

/ 1 n T ' m (S?i n T > m 

V _ At +Vh ' ( >> 
1 C 1 

+ At* n + At^ At h Uq ° }± h\ K 



e n t 



m+l 



(C.2) 



k a e {i,e} 



Finally, we consider the constraints (14.241) to make these 2 equations uncoupled and, 
up to some linear combinations, they can be rewritten under the form 



(v fc - ((V^' m+1 )™ +1 )) + A 



1 rn r ' m+1 | A . = ri? T ' m+1 | A . , 



and 



- (v, • (y> m+i (v 

with Ai, A 2 , iT' m+1 | K and 5 r ' m+1 k defined by 

1 + e 



+ A 2 r0 T ' m+1 k =r5' 
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e-T e 
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